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Abstract. Field theory tools are applied to analytically study fluctuation and 
correlation effects in spatially extended stochastic predator-prey systems. In the mean- 
field rate equation approximation, the classic Lotka- Volterra model is characterized 
by neutral cycles in phase space, describing undamped oscillations for both predator 
and prey populations. In contrast, Monte Carlo simulations for stochastic two-species 
predator-prey reaction systems on regular lattices display complex spatio-temporal 
structures associated with persistent erratic population oscillations. The Doi-Peliti 
path integral representation of the master equation for stochastic particle interaction 
models is utilized to arrive at a field theory action for spatial Lotka- Volterra models 
in the continuum limit. In the species coexistence phase, a perturbation expansion 
with respect to the nonlinear predation rate is employed to demonstrate that spatial 
degrees of freedom and stochastic noise induce instabilities toward structure formation, 
and to compute the fluctuation corrections for the oscillation frequency and diffusion 
coefficient. The drastic downward renormalization of the frequency and the enhanced 
diffusivity are in excellent qualitative agreement with Monte Carlo simulation data. 
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1. Introduction 



In the past two decades or so, analytical and computational tools developed in statistical 
physics have been quite successfully applied to mathematical problems in ecology and 
population dynamics, with the overall goal to arrive at a quantitative understanding 
of the emergence of biodiversity in nature [4]. The typical physics approach to 
complex dynamical systems is of course to first consider perhaps oversimplified idealized 
models that however are designed to hopefully capture the essential phenomenology. 
A considerable part of the mathematical biology literature largely addresses coupled 
deterministic equations of motion for interacting population species that are ultimately 
based on mean-field type rate equation approximations, whereas leaving aside some of 
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the biological complexity provides the opportunity to consistently incorporate stochastic 
fluctuations and spatio-temporal correlations, whose crucial importance has long been 
recognized in the field [5]. 

This paper addresses predator-prey competition models that are defined via 
reaction-diffusion systems on a regular d- dimensional lattice, and whose rate equations in 
the well-mixed mean-field limit reduce to the two coupled ordinary differential equations 
originally introduced independently by Lotka [6J and Volterra [7] nearly a century ago. 
These stochastic spatial predator-prey models have served as paradigmatic examples 
for the emergence of cooperative steady states in the dynamics of two competing 
populations [8] [10] (see also Ref. [H] for a fairly recent overview). The deterministic 
Lotka-Volterra rate equation model is characterized by a neutral cycle in phase space, 
describing regular undamped nonlinear population oscillations with the unrealistic 
feature that both predator and prey population densities invariably return to their initial 
values. In contrast, computer simulations of sufficiently large stochastic Lotka-Volterra 
systems yield long-lived erratic population oscillations p2] [19], whose persistence 
can be understood through a resonant stochastic amplification mechanism [20] that 
drastically extends the transient time interval before any finite system ultimately 
reaches its absorbing stationary state, where the predator population becomes extinct 
|21|. [22] . In spatially extended systems, the mean-field Lotka-Volterra reaction-diffusion 
equations allow for traveling wave solutions [23]-|25j. In the corresponding stochastic 
lattice realizations, these regular wave crests become spreading activity fronts [26] 
that further enhance both populations' life span, and furthermore induce short-ranged 
but significant positive correlations between representatives of either species, and anti- 
correlations between the predator and prey populations [HI [27]. Over the past years, we 
have investigated various different variants of such two-species stochastic spatial Lotka- 
Volterra models for competing predator-prey populations, and found these intriguing 
spatio-temporal structures to be remarkably stable with respect to modifications of 
the detailed microscopic interaction rules [27] [28] . Even in the presence of quenched 
spatial disorder in the reaction rates, the qualitative features of spatial stochastic Lotka- 
Volterra models remain unchanged, although quite remarkably both the predator and 
prey populations benefit from such environmental variability [29J. 

Many qualitative features of stochastic spatial predator-prey systems are adequately 
captured by the associated coupled mean-field rate equations, augmented with diffusive 
spreading. However, one observes strikingly strong quantitative renormalizations of, 
e.g., the characteristic population oscillation frequency, whose numerically determined 
values in various systems were found to be reduced as compared with the (linearized) rate 
equation predictions by factors in the range 2 ... 6, depending on the reaction rates, both 
in the presence and absence of site occupation number restrictions for the predator and 
prey populations [THE?]. I n addition, the neutral cycle oscillations in the original Lotka- 
Volterra rate equations are undamped; in contrast, when a finite carrying capacity for the 
prey species is imposed, the system relaxes to a stable coexistence fixed point. However, 
starting from random initial states, Monte Carlo simulations for the corresponding 
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stochastic lattice models yield initially damped population oscillations in the coexistence 
phase, even in the absence of any restrictions on the site occupation numbers, i.e., for 
infinite local carrying capacities. These are associated with striking spatio-temporal 
structures, namely spreading waves of prey closely followed by predators. In our Monte 
Carlo simulations, we measured the front speed to be markedly enhanced with respect 
to the mean-field prediction [29] . 

The aim of this present paper is to provide a qualitative and even semi-quantitative 
explanation for these intriguing observations. The Doi-Peliti coherent-state path 
integral representation of the master equation for stochastic interacting particle systems 
[3Q]-[3l] (for recent reviews, see Refs. [351 f36j ) . augmented with a means to incorporate 
restricted site occupation numbers [37], will be employed to gain a comprehensive 
understanding of fluctuation and correlation effects in the thermodynamic limit of 
stochastic spatial predator-prey models. More specifically, a perturbative loop expansion 
to first order in the nonlinear predation rate will be constructed; it will allow us to 
demonstrate the instability of the spatial stochastic system against dynamic structure 
formation, and enable the computation of the fluctuation-induced renormalizations of 
the population oscillation frequency and diffusion coefficient [38J. 




This very same formalism was already utilized in Ref. [TT] to demonstrate with 
the aid of renormalization group arguments that the effective critical field theory in the 
vicinity of the predator extinction threshold that emerges at low predation rate for finite 
prey carrying capacity can be mapped onto Reggeon field theory which encapsulates 
the universal scaling behavior of critical directed percolation clusters [39] [12]. Indeed, 
since the predator extinction transition represents a continuous nonequilibrium phase 
transition from an active stationary to an inactive, absorbing state (in the absence of 
any conserved quantities and quenched disorder), one would expect it to be governed by 
the prominent directed percolation universality class [2] [IE]. There exists now ample 
numerical evidence that the critical exponent values at or near the predator extinction 
transition in spatially extended Lotka-Volterra systems are in fact those of directed 
percolation [9]-[T9]. 

The present work also complements and transcends the treatment in Ref. [UJ where 
the same mathematical framework was utilized as a starting point for a van Kampen 
system size expansion, demonstrating on the Gaussian fluctuation level the persistence 
of population oscillations in the species coexistence phase of stochastic lattice Lotka- 
Volterra models, thus generalizing the zero-dimensional analysis in Ref. [20J to spatially 
extended systems. Here, the fluctuation-induced renormalizations of the characteristic 
population oscillation frequency, damping, and diffusivity will be computed to first order 
in a perturbative expansion with respect to the nonlinear predation rate. It should be 
noted that in contrast with the powerful van Kampen system size expansion, the Doi- 
Peliti formalism captures fluctuations and intrinsic reaction noise in the thermodynamic 
limit, and has been successfully applied to systems governed by strong correlations for 
which a simple Kramers-Moyal expansion and Fokker-Planck truncation fails (see, e.g., 
Ref. [36]). The field-theoretic approach has also been employed as an efficient route to 
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construct an expansion about the thermodynamic limit in spatially extended predator- 
prey systems in Refs. [MJ H9] . 

This paper is structured as follows: The following section begins with a concise 
review of the properties of the Lotka-Volterra mean-field rate equations, including 
the modifications induced by a finite prey carrying capacity, and some crucial features 
observed in Monte Carlo simulations for stochastic two-species predator-prey models on 
a regular lattice. Next the construction of the Doi-Peliti path integral representation is 
explained, and its utility demonstrated by a brief summary of the crucial steps that allow 
a mapping of the Lotka-Volterra system with finite local prey carrying capacity near 
the predator extinction threshold onto Reggeon field theory that governs the directed 
percolation universality class. Subsequently, this formalism is employed to construct 
a systematic perturbation expansion with respect to the nonlinear predation rate in 
the species coexistence phase. We then establish the presence of structure formation 
instabilities, and proceed to compute the renormalized population oscillation frequency 
and diffusivity to one-loop order, and compare our results with simulation data. The 
conclusion summarizes these novel results, and gives an outlook to future investigations. 
Two appendices provide additional technical details and an integral table. 

2. Stochastic lattice Lotka-Volterra models 

We begin by first defining the stochastic interacting particle model under consideration 
through a set of coupled irreversible 'chemical' reactions, and then provide a summary 
of its basic features as obtained in the mean-field rate equation approximation. Next we 
discuss the crucial numerical observations from the extensive literature for stochastic 
spatially extended two-species predator-prey systems. 

2.1. Model variants and mean- field description 

We consider a system comprised of two distinct particle species that propagate diffusively 
with continuum diffusion constants Da/b an d undergo the following stochastic reactions: 



The 'predators' A decay or die spontaneously at rate [i > 0, whereas the 'prey' B produce 
offspring with rate o > 0. In the absence of the binary 'predation' interaction with 
rate A, the uncoupled first-order processes would naturally lead to predator extinction 
a(t) = a(0) e _Ati , and Malthusian prey population explosion b(t) = b(0) e at ; here a(t) 
and b(t) respectively indicate the A(B) concentrations or population densities. The 
predation reaction constitutes a nonlinear interaction that simultaneously controls the 
prey particle number and allows the predators to multiply, thus opening the possibility 
of species coexistence through competition. 



A — > with rate /i, 

A + B -> A + A with rate A' 



(1) 



B -> B + B 



with rate a. 
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In the simplest spatial realization of this stochastic reaction-diffusion model, 
both particle species are represented by unbiased random walkers on a c/-dimensional 
hypercubic lattice (with lattice constant do), and one allows an arbitrary number of 
particles per lattice site (see Ref. |27J). The reactions (JT|) can then all be implemented 
strictly on-site: Offspring particles are placed on the same lattice point as their parents, 
and the predation reaction happens only if an A and a B particle meet on the 
same lattice site. If one further assumes the populations to remain well mixed, and 
consequently ignores both spatial fluctuations and correlations, the coupled reactions 
(JTJ) can approximately be described through the associated mean-field rate equations 
for spatially homogeneous concentrations ait) = (a(x,t)), b(t) = (b(x,t)), where a(x,i) 
and b(x, t) respectively denote the local predator and prey densities. With A = a$ A', 
this leads to the two classical Lotka-Volterra coupled ordinary nonlinear differential 
equations jl]: 

a(t) = Xa(t)b(t) - fia(t) , b(t) — a b{t) — A a(t) b(t) . (2) 

The rate equations ([2} display three stationary states (a s ,b s ): (i) the empty 
absorbing state (total population extinction) (0, 0), which is obviously linearly unstable 
if cr > 0; (ii) an absorbing state wherein the predators go extinct and the prey population 
diverges (0, oo), which for A > is also linearly unstable; and (iii) a species coexistence 
state (a u = a/A, b u = /i/A), which represents a marginally stable fixed point with purely 
imaginary eigenvalues ±iuo of the associated Jacobian stability matrix, with the (linear) 
oscillation frequency uq = In /o = y/JHo about the center fixed point (a u ,b u ). In the 
full nonlinear ordinary differential equation system (jSJ), the phase space trajectories 
are determined by da/db = [a (A b — n)]/[b(a — A a)], for which one easily identifies 
the conserved first integral K{t) = X[a(t) + b(t)] — oTna(t) — /ilnfr(t) = K(0). As a 
consequence, the solutions of the deterministic Lotka-Volterra rate equations form closed 
orbits in phase space that describe regular periodic nonlinear population oscillations 
whose amplitudes and phases are fixed by the initial configuration. Naturally, the 
neutral cycles of the coupled mean-field rate equations (j2J) that appear independent of 
the set rates and take the system precisely back to its initial configuration after one 
period represent biologically unrealistic features, and are moreover indicative of the 
fundamental instability of this deterministic mathematical model with respect to even 
slight modifications jl]. 

One important example of such an alteration that aims at rendering the Lotka- 
Volterra system more relevant to ecology is to introduce a finite carrying capacity 
(maximum total particle density) p > that limits the prey population growth [4J. It 
can be interpreted as originating from, e.g., limited food resources for the prey species. 
Within the mean-field rate equation framework, the second differential equation in (j2j) 
then becomes replaced with 

b{t) = a b(t) [1 - b(t)/p) - A a{t) b{t) . (3) 

Again, one finds three stationary states in this restricted Lotka-Volterra model: (i) total 
extinction (0,0); (ii') predator extinction and prey saturation (0,p), which becomes 
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linearly stable for sufficiently small predation rates A < A c = p/p; (iii') species 
coexistence (a r , b r ) with a r — (1 — p/ pA) cr/A and b r = p/A, which comes into existence 
and is linearly stable if the predation rate exceeds the threshold A c . The finite carrying 
capacity causes the Jacobian matrix eigenvalues to acquire negative real parts, indicating 
an exponential approach to the stable fixed point (a r ,b r ): 



p a 



, ± ,,i_l£*f£*_, 



0~ \ p 



(4) 



The neutral cycles of the unrestricted model (j2J) are thus replaced either by a stable 
node for which e± are both real, namely when a > o~ s = 4Ap(pA/p — 1) > 0, or 
alternatively p/p < A < A s = (1 + yl + a/p) p/2p; or by a stable focus with complex 
conjugate stability matrix eigenvalue pair, and consequent spiraling relaxation towards 
the fixed point if a < a s or A > X s . In this situation, both predator and prey populations 
approach their stationary values (a r , b r ) via damped oscillations. Adding spatial degrees 
of freedom, finite local carrying capacities can be implemented in a lattice model through 
limiting the maximum occupation number per site for each species. Most drastically, 
one can permit at most a single particle per lattice site (as, for example, implemented in 
Ref. [H]); the binary predation reaction then has to occur between predators and prey 
on nearest-neighbor sites, and new offspring is to be placed on adjacent positions. In that 
case, one may actually dispense with hopping processes, since all particle production 
reactions automatically entail population spreading. Upon adding diffusive spreading 
terms (with diffusivities D^, Db) to the mean-field rate equations, one may describe 
spreading activity fronts of prey invading empty regions followed by predators feeding 
on them. A well-established lower bound for the front propagation speed is [23j EH H] 



^front > V 4 D A (A p - p) . (5) 

To summarize, within the mean-field rate approximation, a finite prey carrying 
capacity p, which can be viewed as the mean result of local restrictions on the prey 
density originating from limited resources, crucially changes the phase diagram: There 
emerges an extinction threshold (at A c for fixed p) for the predator population, which 
in a spatially extended system in the thermodynamic and infinite-time limit becomes a 
genuine continuous nonequilibrium phase transition from an active to an absorbing state. 
In addition, deep in the species coexistence phase the restricted Lotka-Volterra model is 
characterized by transient decaying population oscillations, which become overdamped 
upon approaching the predator extinction threshold. 



2.2. Monte Carlo simulation results in the species coexistence phase 

Various authors have studied individual-based stochastic lattice predator-prey models, 
predominantly in two dimensions and typically with periodic boundary conditions, that 
in the well-mixed mean-field limit reduce to the classical Lotka-Volterra system; see 
Refs. [8] [19] and [27] for a partial listing. The following is a concise summary of some 
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fundamental results from these extensive numerical investigations, as pertinent for the 
subsequent field-theoretic analysis. 

Monte Carlo simulations in two dimensions, both in the absence and presence of 
local density limitations, observe the emergence of prominent spatio-temporal structures 
associated with remarkably strong fluctuations in the species coexistence phase, even 
far away from the continuous nonequilibrium predator extinction transition. Spherically 
expanding growth fronts of prey closely followed by predators periodically sweep the 
system; any small surviving clusters of prey then serve as nucleation centers for 
new population waves that subsequently interact and for large population densities 
eventually merge with each other. These spreading activity fronts are especially sharp for 
the site-restricted model variants, whereas in simulation runs performed with arbitrarily 
many particles per site, the fronts appear more diffuse and localized [26]. Equal-time 
density correlation functions can be employed to determine the spatial width ~ 10 ... 20 
lattice sites of the spreading activity regions [TTJ [27J [29] . In comparison with the mean- 
field bound (jSJ), the front velocity was measured to be typically enhanced by a factor 
up to ~ 2 ... 3 in simulation runs starting with a single localized activity seed [29J . 

Averaging over the weakly coupled and periodically emerging structures yields long- 
lived but damped population oscillations. As the system size increases, one observes 
the relative oscillation amplitudes to decrease; in the thermodynamic limit, the quasi- 
periodic population fluctuations eventually terminate. Yet locally density oscillations 
persist for both predators and prey species. In the absence of spatial degrees of freedom, 
these can be understood by performing a van Kampen system size expansion about 
the absorbing steady state [20]. The fluctuation corrections may then essentially be 
described by means of a damped harmonic oscillator driven by white noise that will on 
occasion resonantly incite large-amplitude excursions away from the stable fixed point 
in the phase plane. 

From the prominent peaks detected in the Fourier-transformed concentration 
signals, characteristic oscillation frequencies can be inferred [TTJ [271 29J • The thus 
numerically determined typical population oscillation frequencies are found to be 
reduced by a factor ranging between 2 and 6 (depending on the other rates) in 
the stochastic spatially extended system as compared to the mean-field prediction, a 
considerable downward renormalization obviously caused by fluctuations and reaction- 
induced spatio-temporal correlations; compare Fig. 9 in Ref. [11] and Fig. 6(b) in 
Ref. [27]. However, the measured oscillation frequencies / roughly follow the square-root 
dependence on the rates /i and a suggested by the linearized mean-field approximation: 
ujq = 27t/ = y/Jj^, yet with noticeable deviations once either a or /x significantly 
differ from unity. In addition, the functional dependence of / on the rates /i and a is 
surprisingly similar, at least in a mid-range interval of values for both rates near 1. As 
we shall see in section HJ these observations and quantitative trends are remarkably 
accurately reproduced by a first-order analytic perturbation theory for fluctuation 
corrections in stochastic lattice Lotka-Volterra models. 

As the predation efficiency A is reduced (with all other parameters held constant), 
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stochastic lattice Lotka-Volterra systems with site occupation restrictions display just 
the same qualitative scenarios as revealed by the mean-field analysis for eqs. (j3]) with 
finite prey carrying capacity: First, the focal stationary points in the phase plane 
are replaced by stable nodes (corresponding to real stability matrix eigenvalues); the 
population oscillations then cease, and no interesting spatial structures are discernible 
aside from hardly fluctuating localized clusters of predators in a 'sea' of prey that almost 
fill the entire lattice [26]. At a sufficiently small critical value A c , at last the predator 
extinction threshold is encountered, and the measured critical scaling laws near this 
active- to absorbing state transition are very well described by the accepted critical 
exponents of directed percolation [9]-[T9]. 

Simulations in one dimension (usually on a circular domain) yield a crucial 
difference between model variants that incorporate or neglect site occupation number 
restrictions: In the former situation, the A and B particles quickly segregate into distinct 
domains, with the predation reactions occurring only at their boundaries. The long- 
time evolution is consequently dictated by the very slow coarsening of merging predator 
domains [TTj . In contrast, in the absence of site occupation restrictions, one observes 
the system to invariably remain in an active fluctuating coexistence state [27] . 

We finally remark that the above statements naturally all pertain to sufficiently 
large lattices. Of course, any finite system with an absorbing steady state will in 
principle eventually reach and remain trapped in it. However, the associated tpyical 
extinction times are understood to grow fast with system size, namely according to a 
power law [211 E2] ; simulation runs performed in reasonably large lattices consequently 
never reach this extinction state during their entire duration. 

3. Field-theoretic analysis 

This section will first provide a brief overview how a coherent-state path integral 
representation can be constructed directly from the fundamental master equation that 
defines a stochastic interacting particle system [30] [36], see also Refs. [501 [51]; R- e f- [M] 
provides an illustrative alternative derivation. This field-theoretic representation 
faithfully encodes statistical fluctuations, including those caused by discreteness and 
the internal reaction noise, as well as emerging correlations in spatial reaction-diffusion 
systems, and allows for systematic approximative analysis, as will be detailed below 
for two-species predator-prey models. For the sake of completeness, the essential 
steps of mapping spatial stochastic Lotka-Volterra models with site occupation number 
restrictions near the predator extinction threshold onto Reggeon field theory [11] 
will be repeated here as well. The subsequent section H] is then concerned with 
fluctuation corrections in the two-species coexistence phase, which become manifest 
through propagator renormalizations. 
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3.1. Field theory representation 

The Doi-Peliti approach is based on the fact that at any time the configurations in 

locally reacting particle systems can be enumerated through specifying the occupation 

numbers of each species per lattice site i, say here rij for the predators A and m$ for the 

prey B, and that the effect of any allowed stochastic process is to merely modify these 

on-site integer occupation numbers. For now, arbitrarily many particles of either species 

are allowed to occupy any lattice point: rii, rrii = 0,1, ... ,00. The master equation for 

our local reaction scheme ([1]) that governs the time evolution of the configurational 

probability to find n, predators and m t prey on site % at time t through the balance of 

gain and loss terms reads 

d r 1 

— P(n h mi, t) = fj, (m + 1) P(ui + 1, m^, t) - n { P(n h m,; t) 
at - J 

+ 0- 



mi - 1) P(m, mi - I; t) - mi P(m, m,; t) (6) 
- 1) {mi + 1) P(ni - I, mi + I; t) - m P(m, m-t) . 



As initial condition, we may for instance choose a Poisson distribution P(ni, mi] 0) = 
n T Q m^ 1 * e~ no ~ m ° /nil with mean initial predator and prey concentrations n and m . 

Because all reactions just change the site occupation numbers by integer values, a 
Fock space representation is particularly useful. To this end, we introduce the bosonic 
ladder operator algebra [oj, %] = 0, [a i; 0]] = Sij for species A, from which we construct 
the predator particle number eigenstates |nj), Oj jn,) = \rii — 1), a\ \m) = |rij + 1), 
a\ai\ni) = rii |n«). A Fock state with n« particles on site 2 is obtained from the 
empty 'vacuum' configuration |0), defined via Oj |0) = 0, through \m) = ap\0). In 
the same manner, we proceed for the prey particles, with associated annihilation and 
creation operators \ and b\ that all commute with the predator ladder operators: 



\a 



t ,bj] = = [ai,b]]. 



To implement the stochastic kinetics for the entire lattice, one considers the master 
equation for the configurational probability P({ui}, {mi}; t), given by a sum over all 
lattice points of the right-hand side of eq. (jSJ), and recognizes that a general Fock state 
is constructed by the tensor product |{rij}, {mj}} = FT,- |rij) \rrii). One then defines a 
time-dependent formal state vector through a linear combination of all possible Fock 
states, weighted by their configurational probability at time t: 

|$(f)>= Yl P({ni},{m};t)\{ni},{mi}) . (7) 

This superposition state thus encodes the stochastic temporal evolution. Straightfor- 
ward manipulations now transform the time dependence from the linear master equation 
into an 'imaginary-time' Schrodinger equation, governed by a time- independent stochas- 
tic Liouville time evolution operator H: 

9^p-=-H\*{t)), or |*(f)> = e- Ht |*(0)> . (8) 
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For on-site reactions, H ieac = £\ H r (al, &}; Oj, bi) is a sum of local (normal-ordered) 
contributions; for the Lotka-Volterra predator-prey system one obtains 

#reac = - ^ (1 - a}) a f + a (&J - 1) ftj h + A' (aj - b\) a\ a, ft*] . (9) 

i 

Note that each reaction process is represented by two contributions, originating 
respectively from the gain and loss terms in the master equation. For nearest-neighbor 
hopping of particles A(B) with rate D' A (D' B ) between neighboring lattice sites (ij), one 
finds the additional contributions 

tfdifr = J2 [ D 'a - «}) («* - «i) + ^ (&I - 6}) & - bj)] • (10) 

<*i> 

Our goal is to compute averages and correlation functions with respect to the 
configurational probability P({nj}, {m^}; t) which is accomplished by means of the 
projection state (V\ = <0| Ili e ai e b % for which (P\Q) = 1 and (V\a\ = (V\ = (V\b\, 
since [e a %aj] = e ai Sij. For the desired statistical averages of observables O, which 
naturally must all be expressible as functions of the occupation numbers ra» and m*, one 
obtains 

(0{t))= 0({n l },{m t })P({n t },{m t };t) = (V\0({a\a l },{b\b t })mt)) ■ (H) 

{ni},{m;} 

As a consequence of probability conservation, one finds for = 1: 1 = (P|$(t)) = 
( , P|e _ - H ' t |$(0)). Therefore (7- > |i? = must hold; upon commuting e^^ ai+6i ^ with if, the 
creation operators are effectively shifted by 1: a] — > 1 + aj, fe| — > 1 + bj. The probability 
conservation condition is thus satisfied provided Hi(a\ — > 1, fet — > 1;^,^) = 0, which 
is of course true for our explicit expressions © and (IIP]) . Through this prescription, 
we may replace aja, — >■ a« and t»J 6j — >■ 6, in all averages; e.g., the predator and prey 
densities become a(t) = (ctj(t)) and b(t) = (bi(t)). 

In the bosonic operator representation above, we have assumed that no restrictions 
apply to the particle occupation numbers rij on each site. If < 2s + 1, one may 
instead employ a representation in terms of spin s operators. An alternative approach, 
devised by van Wijland, utilizes the bosonic theory, but incorporates site occupation 
restrictions through explicit constraints, which ultimately appear as exponentials in the 
number operators [37]. For example, limiting the local prey occupation numbers to 
or 1 modifies the birth process in fl9j to H ia = a (1 - &t) b\ bi e" 6 « b \ Instead, one could 
also just add a reaction that restricts the local population numbers, e.g., B + B — > B 
with rate z/, yielding an additional term H iu i = —v' (1 — b\) b\ bf. 

As a next step, we follow a well-established route in quantum many-particle theory 
[52] and proceed towards a field theory representation via constructing the path integral 
equivalent to the 'Schrodinger' dynamics (JHJ) based on coherent states, which are defined 
as right eigenstates of the annihilation operators, Oj |atj) = ai \a>i) and bi = (3i |/3j), 
labeled by their complex eigenvalues and fa. One readily confirms the explicit formula 
|aj) = exp(— \ |«i| 2 + ajo|)|0), the overlap integral {aAa^j = exp(— \\ai\ 2 — \\otj\ 2 + 
a*ai), and the (over-) completeness relation J FJ i d 2 a.i = ir. Splitting the 
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temporal evolution (JSJ) into infinitesimal increments, inserting the (over-)completeness 
relation at each time step, and further straightforward manipulations (details can be 
found in Ref. [36]) eventually yield an expression for the configurational average 

(0(t))<x j\[da i da*df3 i df3*0({a i },{f3 i }) exp(-S[a*, a u ft; t]) , (12) 
with an exponential statistical weight that is determined by the 'action' 

S[a*, 0*; a i: ft; t] = £ j^t' (a* ^ + ft ^ + # P (a£, /?*; a,, ft) 



Oi(*) - ft(t) - n a*(0) - m ft*(0) 



(13) 



where the second term at the final time t stems from the projection states, while the 
last one originates in the initial Poisson distributions. Through this procedure, in the 
original quasi-Hamiltonian the creation and annihilation operators a\(b\) and cti(bi) are 
at each time instant replaced with the complex numbers a* (ft*) and aij(ft). 

Finally, we proceed to take the continuum limit, Y2i ~~ ^ a o d I d d x, onit) — > a d a(x, t), 
(3i(t) — > Oq b(x, t), where a denotes the original microscopic lattice constant, whereupon 
the continuous fields a and b acquire dimensions of particle densities, and a*(t) — > d(x, t), 
(3*(t) — > b(x,t), such that a and b remain dimensionless. The 'bulk' part of the action 
then becomes 

° - D A V 2 ) a + b(^- D B V 2 ) b + H r (a, b; a, b) 



S[d } b; a, b] = Id x dt 



(14) 



dt J \dt 

where the discrete hopping contribution ffTUl) has turned into a continuum diffusion 
term, with diffusivities Da/b = a oD' A ^ B . We have thus arrived at a (mesoscopic) 
field theory for stochastic reaction-diffusion processes, with its dynamics governed by 
two independent fields for each particle species, without invoking any assumptions on 
the form of the internal reaction noise. For the Lotka-Volterra reactions ([T]) with 
site occupation number restrictions and/or population-limiting reactions with diffusive 
spreading in d spatial dimensions, the bulk action ( 1T4"|) reads explicitly [UJ 



S[a, b;a,b]= d x dt 



dt 



-D A V 



ot 



D R V 2 6 



+ fi(a-l)a-a(b-l)bbe- a ° bb + v(b-l)bb 2 -\(a-b)aab 

with v — £Jq v' and A = Oq A'; for unrestricted site occupation numbers, the exponential 
term just needs to be replaced with 1, and v set to 0. Expanding e~ a ° bb ~ 1 — a^bb in 
the limit a$ — > effectively replaces the 'hard' exponential constraint with a 'softened' 
particle number restriction, which will henceforth be used. The action (I15j) may now 
serve as a basis for further systematic coarse-graining, constructing a perturbation 
expansion as described below, or, if required, a subsequent renormalization group 
analysis [33 EDI EE] • 

The associated classical field equations follow from the stationarity conditions 
SS/5a = = SS/6b, which are always solved by a = 1 = b, reflecting probability 
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conservation, and SS/5a(x,t) = = 5S/5b(x,t), which yields precisely the mean-field 
rate equations augmented by diffusion terms. Setting a = 1 = b, for a® = = v one 
indeed arrives at the Lotka-Volterra rate equations ([2]), without any restrictions on the 
prey population density, plus diffusive spreading. The modified prey density equation 
03]) with diffusion follows instead, if either v = and a 'soft' particle number restriction 
is implemented with the natural identification p = <2q d , or alternatively with a = but 
adding a pair coagulation reaction with rate v = oj p. It is thus convenient to perform 
the field shift d(x,t) = 1 + d(x,t), b(x,t) = 1 + b(x,t), whereupon the action becomes, 



S[a, b;a,b]= d x dt 



- a b 2 b + - (1 + b) a b b 2 - X (1 + d) (a - b) a b 
P 



, (16) 



with integer a = 2 parametrizing a softened restricted prey occupation, whereas a = 1 
captures instead the presence of the binary reaction B + B — >■ B; the unrestricted model 
is of course recovered for p — > oo. 

We remark that for a — 1, the action ( Tl6l) is equivalent to the two coupled Langevin 
stochastic equations of motion 
da(x, t) 



dt 



(D A V 2 - fi) a(x, t) + A a(x, t) b(x, t) + C(£, t) , (17) 



d 6 ^'*) = m B v 2 + a) b(x, t)-- b(x, t) 2 - X a(x, t) b(x, t) + ri(x, t) , 
at p 

i.e., the diffusive rate equations for the local particle densities, with added Gaussian 

white noise with vanishing means, (0 = 0= (77), and the (cross-)correlations 

(C(x, t) C(x', If)) = 2X a(x, t) b(x, t) 8{x - x') 5{t - f) , 

(C(x, t) r](x', t')) = -X a(x, t) b(x, t) 8{x - x') 5(t - t') , (18) 

(r](x,t)r](x',t')) = 2crb(x,t) l-b(x,t)/p 5(x - x') 5{t - t') , 

describing multiplicative noise terms that vanish with the particle densities, as 
appropriate for the absorbing state at a = = b. Similar Langevin equations were 
derived in Ref. [17] . 

The equivalence of eqs. (j!7p and (lisp with the action ([175]) follows immediately 
from the standard Janssen-De Dominicis field theory representation of Langevin 
dynamics [53] 151] (see also Refs. [50] |5T]). according to which the set of Langevin 
equations dsi(x,t)/dt = Fi[{si(x,t)}] + Q(x,t) with (Q) = and noise correlations 
{d(x, t) Ci(^) t')) = 2Lij[{si(x, t)}} S(x — xf) 5(t — t') is governed by the action 



sm-, { Sl }} = Jd d xjdtj2 - msi}\) - Y^SiLiAM] § 



(19) 



For a = 2, the action ([16]) contains a cubic term of the 'auxiliary' field b, and a 
direct Langevin representation is not obviously possible. In the following, the field 
theory action ffl6]) will serve as the starting point for further manipulations (i) to 
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briefly recapitulate the identification of critical directed percolation as the universality 
class that governs the continuous active to absorbing state transition at the predator 
extinction threshold [IT] , and (ii) to compute the fluctuation- induced renormalization to 
lowest order in the predation rate for the population oscillation frequency and damping, 
as well as the diffusion coefficient in the two-species coexistence phase [38J. 



3.2. Predator extinction transition and Reg g eon field theory 

Here we provide the basic steps by which the effective field theory that describes 
the universal scaling behavior near the predator extinction threshold is constructed, 
following Ref. [UJ. For A ~ A c = p/p, very few predators remain, while the prey almost 
fill the entire lattice, a(x,t) w a r = 0, b(x,t) w b r = p. The reaction scheme ([lj is 
thus essentially replaced with A — > and A — > A + A. We then also need to add 
a growth-limiting process for the predator population, for example again through the 
binary coagulation reaction A + A — > A, say with rate r; heuristically, we have then 
already arrived at the standard single-species death-birth-annihilation reactions that in 
essence define directed percolation processes (see, e.g., Refs. [36l H2| [50], [5Tj ) . 

In the Doi-Peliti representation, we consequently transform the action ( fT6i) to 
new fluctuating prey fields e(x,t) = p — b(x,t) with vanishing mean (e) = 0, and 
e(x,t) = —b(x,t). With the additional predator pair coagulation reaction, this yields 

"~ ' p — A p ] a + t a (1 + a) a 2 



S[a, e; a, e] 




d d x I dt 



+ e 



dt 



-D B V 



a 



dt 



a 



D A V 



'1-eT-l 



a 



e(p-2e) (1-e) 



e) a ee 2 



\p(d 2 + 



1 + a) e ) a + A (1 + a) (a + e) a e 



(20) 



Next we note that the birth rate is a relevant parameter in the renormalization group 

sense, which scales to infinity under scale transformations; this observation simply 

expresses the fact that fluctuations of the nearly uniform prey population become 

strongly suppressed through the 'mass' term oc a for the e fields. It is therefore 

appropriate to introduce rescaled fields <fi(x,t) = y/ae(x,t) and <p(x,t) = t/ae(x,t), 

and subsequently take the limit o — > oo, which yields the much reduced effective action 

d „ _ \ 2 

/i — A p \ a — \ pa a 



Soo[a, 0; a, . 



d a x / dt 



dt 



-D A V z 



rail + a) a + 



+ a pi 



(21) 



Since the fields <fi and <fi only appear as a bilinear form in the action (I2ip . they can 
readily be integrated out, leaving 



SLIM] 



d d x / dt 



( d 



^j\— + D A [r A -V 



Ulft (if) — ?p ) ?p + T ?p 2 -ip 2 



(22) 



where ip(x,t) = a(x, t) a/t/A p, ip(x,t) = a(x, t) ^JXp/r, r A 
u = y/r\p. This new effective nonlinear coupling u becomes dimensionless at d. 



(p — X p)/Da, and 
4, 



Population oscillations in Lotka-Volterra models 



14 



signifying the upper critical dimension for this field theory. Near four dimensions, 
the quartic term oc r constitutes an irrelevant contribution in the renormalization 
group sense and may be omitted for the analysis of universal asymptotic power 
laws at the phase transition. The action (1221) then becomes identical to Reggeon 
field theory, which is known to describe the critical scaling exponents for directed 
percolation [39] |42]. This mapping to Reggeon field theory hence confirms the 
general expectation that the predator extinction threshold is governed by the directed 
percolation universality class [HI [10] , p3] PS], [T8], [19] , which features quite prominently 
in phase transitions to absorbing states |4TJ H3], even in multi-species systems |45j. 
The universal scaling properties of critical directed percolation are well-understood 
and quantitatively characterized to remarkable accuracy, both numerically through 
extensive Monte Carlo simulations and analytically by means of renormalization group 
calculations (for overviews, see Refs. [441 141)1 142] ). 



4. Fluctuation corrections in the coexistence phase 

We now proceed to investigate and analyze the effect of intrinsic stochastic fluctuations 
and spatial correlations in the two-species coexistence phase (and in the thermodynamic 
limit), by means of a systematic perturbation expansion about the (undamped) mean- 
field theory with infinite prey carrying capacity, p — > oo. Various additional technical 
details are deferred to the three Appendices. 



4-1. Doi-Peliti action in the two-species coexistence phase 



In order to address fluctuation corrections in the predator-prey coexistence phase |38j , we 
start again from the Doi-Peliti field theory action ( !T6|) . and introduce proper fluctuating 
fields c(x, t) = a(x, t) — (a) and d(x, t) = b(x, t) — (b) with vanishing mean: 

a(x,t) = + + c(x,t) , b(x,t) = ^(l + B c ) + d(x,t) . (23) 

A V p A / A 

Here, the mean-field values for the stationary densities have been taken into account 
already, such that the counter-terms A c and B c , which are naturally determined by 
the conditions (c) = = (d), contain only fluctuation contributions; this is in accord 
with standard procedures for perturbation expansions in ordered phases [55] [57]. Upon 
inserting (1231) into ()16p . and renaming a(x,t) = c(x,t) and b(x,t) = d(x,t), one arrives 
at the action S[c, d; c, d] in terms of the new fields. It is a sum of three contributions, 



S s [c, d; c, d] 



a p 



d d x / dt 



Bj 1 



_p_ 

pX 



+ A 



+ a + B c 



p_ 

pX 

/' 



[1 + B C ) [A c + ^B c )d+[l-^- + A c ) (l + B c )c(c-d) 



l-a-z- 1 

px 



B r 



d 2 



[a 



pX 



1) ^- (1 + B c f d 3 



(24) 
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which represent source terms, the bilinear or harmonic contributions 



Sh[c, d; c,d]= / d d x I dt 



pX 



r[ ^--D A V 2 -pB c )r (25) 



a [I — ^ + A c ) cd + d [ — D B V 2 



and finally the nonlinear vertices 



S v [c, d; c, d) — — J d d x J dt 



d_ 

dt 



p(l + B c 



a 



A c +-^-(l + 2B C ) 
pX 



d 



c(c 



l> 



+ a (1 - + A c ) c(c- d)d + a 
pX 



d) c 



l> 



l-2a^-{l + B c 
pX 



d 2 d 



(26) 



- 2(a - 1) ^ (1 + B c ) d 3 d + X (1 + c) (5 - d) cd - - (1 + d) a dd 2 
pX p 

(recall that the exclusion parameter assumes only the values a = 1 or 2). Note that 
since the definitions (123 j) already contain the mean-field expectation values of the field, 
the linear source terms ~ c, d in fl24|) are mere counter-terms. 

The integrand in the harmonic action ( l25|) can be written as a bilinear form 
(c c?) A . Defining the spatial and temporal Fourier transform for an arbitrary field 
via 



<f>(x , t) 



d d q 
(2vr) c 



^0(g»e^—) 



(27) 



(and omitting the fluctuation corrections ~ A c , B c ), we have in Fourier space 

A(,.^)= j - . . / , i . (28) 



-iiO + D B q + a fi/p X 



The next step is to diagonalize the non-symmetric bilinear coupling matrix A = A(0, 0). 
To this end, we need its right and left eigenvectors, Ae± = X±e±, f± A = \± f± 
that satisfy the orthogonality relation /J e T = 0. Introducing the eigenvector matrices 
P = (e+ e_) and Q = (/+ /_), one then readily confirms Q T AP = diag(A + A_), 
with the diagonal elements A± = X±f±e±. Upon defining new fields (p± and (p± via 
CD = P (£) ^nd (5 d) = (£+ £_) Q T , finally (c d) I Q = (ft. ft) diag(A + , A_) (£) = 
A_)_(^ + (/9 + + A_(^_(/?_ . The eigenvalues of the matrix A are just the negative of the stability 
matrix eigenvalues in the coexistence phase, A± = ±iaj +7o — — e ±> c -f- ec i- ®; with the 
mean-field ('bare') oscillation frequency uq and damping constant 70 (see also Ref. [17]): 



pA 



7o 2 , 



7o 



p cr 

277a 



(29) 



Observe that u 2 = pa and 70 — > as the carrying capacity p — > 00: There is no damping 
of the mean-field oscillations in the absence of local particle number restrictions; in 
this situation, damping terms are in fact generated by stochastic fluctuations, as will 
be demonstrated below. Choosing the eigenvectors e± = (iuo =F 70 ± p) /iujQt/2/2, 
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f± = (±/i iu Q ± 7o)/iwov / 27t, with f±e± = ±l/iu , the harmonic action (125^) is 
diagonalized by means of the linear field transformations 

1 / 



2/i 
/i £ 



V?+ + <P- ~ 7o ■ 



IU 

1 



d 



2 zwo 
Indeed, upon inserting fl3"0"|) into 



V9+ + + 7o 



2 zcj 



(30) 



2yU V ' 2Co> 

one obtains the harmonic action in terms of 



the new fields 

Sh[<p±;<p± 



+ 



1 

+ 7o ~ ^ 
2zwo 
+ 7o 



d d £ /dt 



aA c + 



- 7o)(iwo + 7o - A 4 ) + 47o(ic;o + 7o 



+ ^^V 2 + ^ + 7o 



iu - 7o 



tu 



D' V 2 



iujQ + 7q - /i 

2zwo 
2W - 7o + A* 



2iu;, 



a Ac 



aA c 



2icu 

(zw + 7o)(iwo-37o-A4) 



(iuj - 7o)(iwo + 37o + //) 



o 



<9 



-Dn V 2 — — -Dn V 



9t 



iw 

(iw + 7o)(^o 



iu + 70 + 



7o + A 4 ) 



2icoo 
iu - 7q + n 
2iu 

47o(iw - 



5, 



B, 



o-A c 



(31) 



7oJ 



2zwo 



5, 



where D = (Da + D b )/2 denotes the mean particle diffusivity, and D' = (Da — D b )/2 
indicates the asymmetry in the diffusion coefficients. In the following, we shall restrict 
ourselves to the case of equal diffusivities Da = D b = D and D' Q = 0; the harmonic 
propagators in the diagonalized theory then read in Fourier space 

±icjQ 



((p±(q,u)(p±((?,u'))o 



(2vr) d+1 5{q + ?) 5(to + to') , (32) 



— iu + D q 2 ± iujQ + 70 

whereas the off-diagonal two-point correlation functions (<p±(q, to) <£>q=(<f, to')) contain 
only counter-terms and hence vanish in the harmonic approximation. Akin to spin waves 
in magnets, the poles of the propagators ( 1321 describe (anti-)clockwise propagating 
waves with frequency u and damping 7 , with additional diffusive relaxation ~ D q 2 . 
The delta functions in (1321) reflect spatial and temporal time translation invariance. 

Upon expressing the sources (j2H) and nonlinear contributions (T261) as functionals of 
the new fields, a multitude of terms is generated which renders any subsequent analysis 
quite cumbersome, see Appendix A Consequently we shall address the limit of large prey 
carrying capacity p — > 00, for which the mean- field approximation predicts undamped 
oscillatory modes with frequency uo = y/JIa, see eq. (jlj). Correspondingly, we shall 
henceforth retain finite p and non-zero damping 70 solely in the propagator terms (I3ip. 
but set = = 70 everywhere else. The source terms then just read 

2 a 



S 3 [(p±](p±] = / d d x dt 



ito A c (l + B c ) -p:B c (l + A c ) 



Population oscillations in Lotka-Volterra models 



17 



+ 



iu A c (l + B c ) +fxB c (l + A c 



l + B r 



+ 2 



2A 



M (1 + A) + o- 



(iu -//)(! + A c ) + cr 



(33) 



(iu + /i) (1 + A c ) - a 



<P- 



compare ( lA.ip in Appendix A The linear source terms are mere counter-terms; following 
eq. f lT9|) . one may interpret the quadratic ones as generated by stochastic noise. 
From eq. f )26|) one obtains the three-point vertices in the limit p — > oo: 



d x dt 



(iu - //) iu (1 +Ac)—fi(l + B c ) 



~2 



(zu; - /i) iw (1 + A c )+/i(l + B c ) 



+ 2U 



2 U 



iw (1 + As) - (1 + S c ) 



iwo (1 + A) + /x (1 + S c 



+ ZCU CT V 2 - 



+ iujq a ip + if- if. 



- (ioj + n) iu (1 +A c )—fi(l + B c ) 



+ + //) iw (1 + A c ) + ii (1 + 5 C ) 



zw a ) (p_(p + 



tu a \ (f?_(p- 



• (34) 



The nonlinear vertices of the full action are listed in eqs. flA.2j) . fl A. 3|) in Appendix A 
Note that in the large carrying capacity approximation, the various field theory 
contributions naturally become independent of the parameter a. Both the reduced and 
full actions remain essentially invariant under exchange of the labels + < — > — , aside 
from complex conjugation, an obvious consequence of the complex conjugate eigenvalue 
pairs A± for A and the corresponding eigenvector symmetry, see eq. (l30l) . Formally, 
this symmetry is conveniently expressed in terms the vertex functions with m± external 
outgoing (p± and n± incoming (p± legs: 



r + m+ (Xj , tj) = T + ™- _m + . + n_ _n + (fj, £j)* 



(35) 



As a direct consequence, T_ 



.m m- in n 



[Si, U) must be real. 



4-2. Counter-terms and propagator renormalization to one-loop order 

The propagators (I3"2"j) along with two two-point noise sources fl3"3"|) and three-point 
vertices (13^|) represent the building blocks for the Feynman diagrams that graphically 
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represent the different contributions in a perturbation expansion with respect to the 
predation rate A, which serves as the nonlinear coupling here [38] . 




Figure 1. Feynman graphs for the expectation values (<p±) up to one-loop order, 
where the '•' symbol in the first diagram represents the counter-terms. 



As our first step, we need to determine the counter-terms A c and B c to first order 
in A from the conditions (<p±) = 0. The associated Feynman graphs are shown in Fig. [H 
and the corresponding analytic expressions read 







+ 



a % 
2 A 



iu Q A c t fJ>B, 
d d k 



\L — <J — IUq 



ji — a + icjo 



A^^jJ J (2ir) d \iu + -fo + D k 2 -iu + 7o + D k 2 



(36) 



These are readily solved, with the result 



A r = B r 



d d k 



l> 



a 



IUq 



l> 



a 



IUJQ 



4lo J (2ir) d \iuj + 7o + D k 2 -iu + 7o + D k 2 



0(A 2 



d d k [i - a + 70 + D k 2 



+ 0(A 2 



(37) 



(2n) d co 2 + ( 7o + D k 2 ) 2 
We may now proceed to the fluctuation renormalization of the propagators (l3"iZ|) to 
first order in the predation rate A. To this end, we require the two-point vertex functions 
r± ; ±((f, oj) to one-loop order. Denoting their fluctuation corrections by r£. ± (<f, u), the 
structure of the low-frequency and small-wavevector expansion is 

7o 

iu iu 



r± ; ±(g,w) 



Re r£? ± (o, o) ± 1°- + i im r« ± (o, 0) ± ± 



dq 2 



ILO 



q=0 



diu 



+ 



u>=0 



Note that the symmetry (I3^|) implies 

r +;+ (g,o) = r_ ; _(g,o)* 



diu 



u=0 



diu 



uj=0 



(38) 



(39) 



The one- loop Feynman diagrams for T±.±(q, u) are depicted in Fig. [2j Performing 
the internal frequency integrals, one arrives at the associated analytic expressions 



r± ; ±(g» 

+ 



±iu 
A (zkiuo — /i) 



iu ± iu + 7o + D Q q 2 



a 



li 



[iul 



[ji{a — 2iuq) =f iu a\ 



± iuoj A c 

(27T)" 



d Jd 



d d k 



iu/2 ± iu + 70 + A) (g 2 /4 + k 2 ) 
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+/- +/- 
Figure 2. Feynman graphs for the vertex functions T±-±(q, lj) up to one-loop order. 
A (±iu - /i) . /■ (2n)- d d d k 



jt/w 2 ' ' '' " 7 iu/2 =f iw + 7o + -Do (<? 2 / 4 + A; 2 ) 

A . w , f d d k 1 

± — ; — [o~ — ji± 2iuo) (a — n± iojQ 



+ 



«*>0 J (27T) d ±ZW() + 7q + £) (g/ 2 + k) 2 

x 1 ±iooo + lo + Do (q 2 /4: + k 2 ) 

±itu + 7o + D Q (q/2 - k) 2 + iu + 7o + D (q 2 /A + k 2 ) 

X , ^ r d d k 1 1 



(v + rfj 



i (2tt)^ 7q + £, (g/ 2 + k) 2 70 + D (q/2 - k) 2 
x 7o + ^o(g 2 /4 + P) 

iw/2T^ + 7o + D (g 2 /4 + A; 2 ) ' V ; 
where the last two terms have been symmetrized with respect to the external wavevector 
q. Naturally, eq. ( HOI satisfies the symmetry constraints (1391) . Clearly, lmr±.±(0, 0) does 
not vanish, which implies that the nonlinear fluctuations generically either generate a 
damping term for the population oscillations, see eq. ( )42|) . or induce an instability 
towards spatial structure formation, as observed in the lattice Monte Carlo simulations. 

Notice furthermore the convolution of both clock- and anti-clockwise propagating 
modes in the 'triangular' fluctuation loop of the last Feynman graph in Fig. |2J As 
a consequence, the imaginary 'mass' terms +%Uq in the first two factors within the 
associated wavevector integral cancel each other, as becomes apparent in the final term 
of eq. ( HUj) . For vanishing damping 7 — > this induces an infrared divergence in d < 2 
dimensions. It is precisely these contributions that cause large fluctuation corrections 
for the renormalized oscillation frequency, eq. ( H6l) below, in the coexistence phase of 
the spatial Lotka-Volterra system even at finite (but small) damping 70 . 

4-3. Renormalized damping, oscillation frequency, and diffusion coefficient 

Appropriate definitions of renormalized oscillation parameters are suggested by the 
functional form ( 1381) of the vertex functions T±-±(q, u). We thus cast the renormalized 
two-point vertex functions in the form 

y\R ( -> \ , , 1R , W Dr q 2 

T±±(q, u) = l±- — ± — ± , (41) 

ZUr Ur IUr 

whence we identify the renormalized damping 7^, frequency lor, and diffusivity Dr via 

7o+^oImrg ± (0,0) 
7fl = : ; TZZ^T) ~ — 7~7~. '■, — » I 42 ) 



l^wolm [dr±! ± (0,u)/diu 



Ld=0 
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D 



R 



g=0 



20 

(43) 

(44) 

lTWolm [dT^ ± (0,oj)/diuj}^ 

Note that a negative 'damping' 7^ < in eq. (141 j) indicates an instability towards a 
spatially inhomogeneous configuration at wavenumber g c = a/|7,r|/-Dr or characteristic 
wavelength A c = 27r^D Q /\^ R \ + 0(A 2 ). 

Upon evaluating the basic one-loop result (HOI , and following the prescriptions 
(Ii2]) - (l44l . it is a straightforward task to compute the renormalized parameters 7^, ur, 
and Dr. Intermediate steps and technical details can be found in Appendix B As a 
final task, one needs to perform the resulting wavevector integrals for the fluctuation 
corrections. With the aid of the integral table in Appendix C one finds with ( 1C.2I) . 
( 1C.7|) . and (1C.9j) for the renormalized or fluctuation-induced damping ( IB.3j) : 



1R 
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= To + A(^) d/2 A7 R + 0(A 2 ). 
The renormalized oscillation frequency ( IB. 51) becomes 
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Notice that the effective expansion parameter in this perturbation series is given by 
(A/wo) (wo/-Do) d//2 ; accordingly we have introduced dimensionless first-order fluctuation 
corrections Atr, Awr, and AD R . Naturally, when diffusion is fast compared to 
the characteristic oscillation time scale, the system becomes well-mixed and spatial 
correlations irrelevant. Deviations from mean- field theory induced by the fluctuation 
loops are then minute. In dimensions d < 2, when we let 70 — >■ 0, the leading 
fluctuation correction to the oscillation frequency diverges ~ (wo/To) 1_d ^ 2 ; it is negative, 
and symmetric under formal rate exchange fx < — > a, c.f. the last term in the 
second line in eq. (1461) . If we interpret 70 in the above equations as a small, self- 
consistently determined damping, these features are in remarkable qualitative agreement 
with our earlier Monte Carlo observations: Fluctuations and correlations induced by the 
stochastic reaction processes induce a strong downward numerical renormalization of the 
oscillation frequency, with very similar functional dependence on the rates fi and a. Note 
that d c = 2 can be viewed as (upper) critical dimension for the appearance of singular 
infrared fluctuation contributions (in the limit of infinite prey carrying capacity p — > 00 
or 70 — » 0), which resemble dynamic coexistence anomalies induced by Goldstone modes 
in systems with broken continuous order parameter symmetry (see, e.g., Ref. [57] and 
references therein). 

In the following, the expressions (145p - (l4"T]) are evaluated at integer dimensions 
d = 1,2,3, and 4. In low dimensions, i.e., for d — 1 and d = 2, the renormalized 
oscillation frequency ( )46|) becomes singular in the limit 70 — > 0, caused by the 
interference of counter-propagating clock- and anti-clockwise internal modes. For the 
renormalized diffusivity D R and the fluctuation- generated damping 7^, these infrared 
singularities cancel out. In d = 1 dimension, the leading terms in 70 are: 
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Figure 3. Fluctuation contributions to the damping Atr, oscillation frequency Awjj, 
and diffusivity AD^ in d = 1 dimension. The fluctuation corrections to the frequency 
depend crucially on the ratio 70 /wo, especially when a C /1 or <r > /i. 
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The dimensionless fluctuation corrections A 7^, AcDr, and AD R , c.f. eqs. (j4"5]) - (j4"7j) . 
are plotted in Fig. [31 The fluctuation-induced contribution to the damping is always 
negative, indicating the instability towards spatially inhomogeneous structures that 
emerge when 70 = A (A'TrI (uq/ Do) d / 2 + 0(\ 2 ). The oscillation frequency is renormalized 
to lower values by the loop corrections, with the leading term ~ a/ w o /to further 
amplified when either o" < /1 or cr > /i. Likewise, fluctuations invariably enhance 
diffusive spreading. The fluctuation corrections all appear remarkably symmetric with 
respect to exchanging a < — > fi, as is evident in Fig. [3] with its logarithmic scale for the 
rate ratio a/fi by the approximate mirror symmetry about the a/fi = 1 axis. 
In d — 2 dimensions, one gets 
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Figure 4. Fluctuation contributions to the damping Atr, oscillation frequency Awjj, 
and diffusivity ADr in d = 2 dimensions. As in one dimension, the fluctuation 
corrections to the frequency strongly depend on the ratio 70 /u>o- 



and the fluctuation contributions are depicted in Fig. HJ The graphs look remarkably 
alike to Fig. [3] for d = 1, but the overall scale of the corrections A^r and ADr is 
reduced by a factor ~ 4, and for Aojr by ~ 10, with its leading term acquiring just a 
logarithmic singularity as 70 — > 0. Again, the system is rendered unstable against spatio- 
temporal structures. According to eq. ()5]), the fluctuation-enhanced diffusivity suggests 
faster front spreading than predicted by the bare mean-field rates, as indeed observed 
in two-dimensional Monte Carlo simulations [29]. The population oscillation frequency 
is strongly renormalized downward, with an approximately equal functional dependence 
on the rates a and fi; moreover, the deviations from the mean-field values grow in 
size as the ratio a/fi is tuned away from unity. These analytic perturbative one-loop 
results are in remarkable qualitative agreement with the Monte Carlo simulation data 
for two-dimensional stochastic Lotka-Volterra systems, as shown in Fig. 9 in Ref. [TT] 
and Fig. 6(b) in Ref. [27]. 

In d = 3 dimensions, one may safely set the bare damping constant to zero, 70 — > 
(or p — > 00) to obtain 
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Figure 5. Fluctuation contributions to the damping A^r, oscillation frequency Awjj, 
and diffusivity A5j{ in d — 3 dimensions. 
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Figure |5] shows the associated fluctuation corrections A^r, Aur, and ADr, which 
compared to one and two dimensions are considerably reduced in magnitude, but 
otherwise display quite similar features. 

In higher dimensions d > 4, the fluctuation corrections become formally ultraviolet- 
divergent, and thus a finite cut-off A in momentum space must be implemented; e.g., in 
d = 4 dimensions: 
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As Fig. [6] demonstrates, the cut-off dependence in the loop corrections is rather weak in 
d = 4 dimensions. For low cut-off values, the fluctuation contributions to the damping 
appear positive in the approximate interval 1 < aj\i < 30, but turn negative in the 
continuum limit of large A, still signaling instability with respect to structure formation. 



Population oscillations in Lotka-Volterra models 



25 



0.00 

,s -°- 04 

< -0.08 

-0.12 
0.00 

-0.04 

'I -0.08 

-0.12 



0.060 
( «l 0.045 
< 0.030 

0.015 

10~ 2 lO" 1 10° 10 1 10 2 

Figure 6. Fluctuation contributions to the damping Ajr, oscillation frequency AQr, 
and difFusivity A.Df> in d = 4 dimensions. Notice the weak logarithmic dependence on 
the ultraviolet cut-off A. 

The typical values of A^, ADr, and Aur are all diminished by factors ~ 4 ... 5 as 
compared to d = 3; the cut-off dependence in the renormalized frequency only becomes 
noticeable for er//x <C 1, see eq. (155)1 . 

5. Conclusion and outlook 

This paper describes in some detail how the stochastic kinetics of spatially extended 
predator-prey systems of the Lotka-Volterra type, as encoded through a classical master 
equation, can be mapped onto a continuum field theory representation, while faithfully 
preserving the internal demographic and reaction noise and the ensuing correlations. 
The connection of the more microscopic Doi-Peliti field theory action with a mesoscopic 
description in terms of coupled Langevin equations was pointed out, and the associated 
white noise correlations were systematically derived. The continuum representation was 
then employed to demonstrate that the predator extinction transition, induced by a 
finite prey carrying capacity, is indeed governed by the universal scaling exponents of 
critical directed percolation, as one would generically expect for such a nonequilibrium 
phase transition from an active to an absorbing state. 

After a brief review of the most striking features of stochastic predator-prey 
models in the species coexistence phase, the Doi-Peliti field theory representation and 
a first-order perturbation expansion with respect to the nonlinear predation rate, in 
the limit of large prey carrying capacity, were employed to qualitatively and semi- 
quant it at ively confirm crucial salient observations from Monte Carlo simulations on 
regular lattices: (i) Spatial predator-prey systems in the species coexistence phase are 
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generically characterized by the emergence of fluctuating spatial structures, namely 
continually expanding and merging activity fronts, (ii) The recurring passages of 
population waves locally incite persistent density oscillations for both predators and 
prey, (iii) Fluctuations in the two-species coexistence phase are remarkably and 
quite unusually strong; as compared with the (linearized) mean-field prediction, they 
considerably renormalize the oscillation frequency, especially in d < 2 dimensions. 
Explicit analytical results for the fluctuation- induced damping, and the renormalized 
oscillation frequency and diffusion coefficient were provided to one-loop order. They 
showed that (iv) the leading fluctuation contribution to the frequency is negative, and 
symmetric in its functional dependence on the rates a and /x; and (v) the diffusivity is 
invariably renormalized upward, implying faster front propagation speeds as compared 
to the mean-field approximation. 

An important open question is which of the numerous standard mathematical 
models in ecology, population dynamics, and chemical kinetics, many of which are 
frequently analyzed merely on the level of mean-field rate equations, are similarly 
strongly affected by stochastic fluctuations and intrinsic correlations. Remarkably, 
and perhaps counter-intuitively, Monte Carlo simulations of stochastic spatial variants 
of cyclic three-species predator-prey systems, namely both spatial rock-paper-scissors 
games (with conserved total population) and the May-Leonard model (which displays 
no global conservation law) do not reveal noticeable fluctuation effects, see Refs. (5HI EH] 
(and further references therein). Apparently, the mechanism causing strong fluctuations 
in the spatial stochastic two-species Lotka-Volterra system, namely the destructive 
interference of counter-propagating internal modes, is conspicuously absent in extensions 
to additional participating species. This fact becomes even more puzzling as the 
stochastic cyclic rock-paper-scissors model has been shown to reduce to the stochastic 
and strongly fluctuating Lotka-Volterra system in a highly asymmetric rate limit where 
a single species becomes abundant [60]. A careful field-theoretic analysis based on the 
Doi-Peliti representation of the corresponding stochastic master equation should be 
capable to shed light on this issue, and hopefully explain this important distinction 
between apparently closely related population dynamics or reaction-diffusion models. 
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Appendix A. Field theory and counter-terms for finite carrying capacity 

In this appendix, we write down the explicit field theory for finite prey carrying capacity 
p, and sketch the evaluation of the associated counter-terms A c and B c to first order in 
the predation rate A. Upon expressing (1241) in terms of the fields (p± and (p± by means 



Population oscillations in Lotka-Volterra models 
of eqs. (1301) . one obtains the source terms 



27 



S s b±; (p±] = jd d xjdt ^ -A ^ (iu + 70) (1 + B c ) (A c + -^B C ^- fiB c x 



1-4 + 4 

pA 

, g(i+g e ; 

2w 2 A 



(iw -7o)(l + 5 c ) A c +-q-B c 1-^r + A 



pA 



(zw + 7o) 2 ( 1 - ^ (1 + S c ) ) - /i(iw + 7o - A*) ( 1 - + A 



pA 



pA 
I' 



pA 



~2 



- 2 



a [i 



pA 



"o + 7o) 1-^(1 + B c ) -p( 7 o-/i) + A 



pA 



+ 



-(«-ih/£ 



(zo; - 7o ) ( 1 - — (1 + -B c ) ) + ^(iwo - 70 + /x) ( 1 



At o- (1 + £ c 



(iw + 7o) + (iwo - 7o) 



pA 

3 



<P- 



(A.l) 



2 2z^pA 2 

Note that the cubic source contributions are absent if a — 1. The nonlinear action (j26l) 
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In addition, there are four- and five-point vertices (the latter arise only for a = 2): 
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However, these do not enter the one-loop analysis, but only contribute to higher orders 
in the perturbation expansion. 
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Figure Al. Feynman graphs for (<p±) up to one-loop order in the full held theory. 

Naturally, with these many contributions in the action, any subsequent perturbative 
calculation becomes quite elaborate and lengthy, as will next be demonstrated by 
computing the counter-terms A c and B c in the full theory. The contributing Feynman 
graphs up to one-loop order are depicted in Fig. IA1I The associated analytic expressions 
for the expectation values (<f±) become to first order in A: 
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Because of the fundamental symmetry (|35|) . separating (1A.4I) into its real and 
imaginary parts yields only two coupled linear equations for A c and B c . By means 
of straightforward (but tedious) algebra one finally obtains 
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In the large prey carrying capacity limit p — > 00 with 70 — > 0, these expression coincide 
and reduce to eq. ([3] 



Appendix B. Evaluation of the one-loop vertex function 

Next we provide some intermediate steps andadditional technical details for the 
evaluation of the propagator self-energy of vertex function T±-±(q,u) that results in 
the renormalized damping coefficient jr, frequency ujr, and diffusivity D R . 

Collecting and rearranging the one- loop contributions in eq. (|40|) . one arrives at 
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It is worth noting that the wavevector integrals are all of order l/k A (or higher inverse 
powers of k) and consequently develop ultraviolet divergences only in dimensions d > 4; 
as they should, the counter-terms have cancelled contributions of order 1/k 2 . From 
eqs. (J4"2l) and flB.2j) one immediately infers the fluctuation-induced damping 
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We furthermore need 
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which along with (IB.lj) provides us with the renormalized oscillation frequency (j43l 
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+ A 
+ A 
+ A 



/i f d d k 



A J (27r) d Wo 2 + ( 7 o + D A ; 2 ) 2 
3(a + /i) 2 /" d d A; 1 



16 



(2ir) d 70 + D A; 2 w 2 + ( 7o + D P) 2 



a 2 + 8 a\i + /i 2 /" d d k 



7o + D A; 2 



16 

[a — /x) a/i f d d k 



16 



(2ny [ w 2 + (7 + J D P)2]2 
1 



(2vr)<* [a; 2 + (7 + DoA; 2 ) 2 P 



(er + /i) 2 (x/i f d d k 



(2ir) d 7o + D A; 2 [a; 2 + ( 7o + D A; 2 ) 
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(a 2 — 4 a/i + /i 2 ) a/i /" c? d A; 



4 



7o + D A; 2 



(2n) d [ w 2 + ( 7o + jDo F) 2 ] 3 



3(a-/i)a 2 /i 2 /" d d fc 



4 



and 
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=F — 1m +— 

^ D dq 2 



q=0 



+A 



+ A 



-A 



(2n) d [ w 2 + (7 + J D P) 2 ]3 
a — fx f d d k 1 



+ 0(A 2 ) , (B.5) 
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(a + /i) 2 f d d k 1 
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8d 
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(27r) d ( 7o + D A; 2 ) 2 u; 2 + ( 7 o + D k 2 f 



-A 



3(a-/i) 2 y d d k 7o + D A; 2 



16 
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^ a 2 - 4 an + /i 2 /" d d A; D k 2 

-A 
+ A 
- A 



8d J (2vr) d [ w 2 + ( 7o + J D A; 2 ) 2 ] 2 
1 1 (a — /i) cr/i /" d d A; 1 



16 J (2n) d [w 2 + (7o + J D A; 2 ) 2 ] 2 
?>{a-n)an f d d k D k 2 ( 7o + D k 2 ) 



2d J (2ix) d [u 2 + ( l0 + D k 2 ) 2 f 
[a + jifaji f d d k 1 1 



16 7 (2n) d 7o + D A; 2 [co 2 + ( 7 o + D k 2 ) 

(a 2 — 4 (x/i + /U 2 ) <T/U ^ d d A; 70 + D k 2 

+ A 
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-A 
+ A 



4 y (2tt)^ [^ + (70 + ^0 A; 2 ) 2 ] 3 

2 - 4 a/i + /i 2 ) (T/i /" d d A; D k 2 



d J (2n) d [ w 2 + (7 + J D A;2)2]3 

3(<T-/i)aV / d d A; 1 

(2^p K 2 + ( 7 o + Dofc 2 ) 2 ] 3 



_ 3 (a - /i) a 2 /i 2 f d d k D k 2 (70 + Dp A: 2 ) 

d J (2rr) d [oo 2 + ( lo + D k 2 )Y ' 



d J (2n) d [ w 2 + (7o + J D P) 2 ] 4 

whence eq. ( )44|) with ( IB.4j) at last yields the renormalized diffusion coefficient 
D R 1 a — /i f d d k 
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Finally, the wavevector integrals for the fluctuation corrections need to be carried 
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Appendix C. Wavevector integrals 



The required integrals are of the following form, and readily evaluated (where convergent 
in the ultraviolet) by means of Euler's Gamma function: 
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(C.4) 



This last result also follows from the sum of eqs. (IC.2p and (IC.3p . if one observes that 
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Next, decomposition into partial fractions gives 
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Taking derivatives with respect to the parameters 70 and/or w one then obtains: 
d d k 70 + -D k 2 1 d f d d k 1 
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For explicit evaluation at d — 2, the Gamma function T(l — d/2) diverges, but its 
poles in the expressions for the renormalized oscillation parameters are all cancelled, 
as can be checked by setting d — 2 — e, and carefully taking the limit e — > 0. Indeed, 
the singularities in two dimensions are eliminated by the counter-terms A c and B c . 
At d — 4, ultraviolet divergences appear, which must be regularized by a cut-off A 
in momentum space that originates from the underlying lattice; e.g., A = 2n/a in a 
hypercubic lattice with lattice constant a . In the above integral listing, these ultraviolet 
singularities emerge as poles in e = 4— d. For example, the logarithmic cutoff dependence 
i ln(l + A 4 DI/uq) is represented in dimensional regularization by T(l + e/2)/e(l — e/2). 
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